%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION IIV_Bounds                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DESCRIPTION:
%   This function implements the inference procedure described on page 666
%   of Nevo and Rosen (2012, Review of Economics and Statistics),
%   henceforth referred to as NR.  It is assumed that the vectors moments_u
%   and moments_l are root n consistent estimators of upper and lower
%   bounds, respectively, on the parameter of interest.
%   The confidence interval produced is a nominal 2p-1 level confidence
%   interval for the population interval bounds based on Bonferroni's
%   inequality, e.g. p=0.975 produces a 0.95 confidence interval for the
%   population interval.
% INPUT PARAMETERS:
%   moments_u: Column vector of upper bound estimates on the parameter of interest
%   moments_l: Column vector of lower bound estimates on the parameter of interest
%   se_u: standard errors of moments_u
%   se_l: standard errors of moments_l
%   Omega_u: sample correlation matrix of upper moments.
%   Omega_l: sample correlation matrix of lower moments.
%       Note: These matrices are decribed as the respective sample variance
%             covariance matrices on NR p. 666, but should be the
%             sample correlation matrices as indicated here.
%       Note: The Matlab chol function gives an *upper* triangular matrix
%             so that if chol(A) = L, then A = L' * L.
%   p: The parameter p_n described on page 666.
%   n: The sample size.
%   c: Constant factor to be used in moment selection if Selection_Method
%       is 'IIV'. The procedure described in NR Section IV is with c=2.
%       If Selection_Method is 'AIS' then adaptive inequality selection is
%       used.
%   R: Number of multivariate normal simulation draws for computing
%       critical values.
%   Selection_Method: Either 'IIV' or 'AIS' to indicate inequality selection
%       based on the procedure on page 666 of NR or
%       moment selection based on adaptive inequality selection as in
%       Chernozhukov, Lee, and Rosen (2013), developed subsequent to NR.
%       The recommended method is AIS, which is the default  if
%       the argument is not 'IIV'.
% OUTPUT PARAMETERS
%   bound_estimates: consistent estimates for the lower and upper bounds
%      on the parameter of interest given by the maximum of the lower
%      bounds and the minimum of the upper bounds.
%   confidence_interval: An nominal 2p-1 level confidence interval for the
%      population interval bounds based on Bonferroni's inequality.
%      For example, p=0.975 produces a 0.95 confidence interval.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bound_estimates,confidence_interval] = IIV_Bounds(moments_u,...
            moments_l,se_u,se_l,Omega_u,Omega_l,p,n,c,R,Selection_Method)

% ENSURE moments_u IS A COLUMN VECTOR
if (size(moments_u,2) > 1)
  if (size(moments_u,1) > 1)
    error('moments_u must be a vector.'); 
  else
    moments_u = moments_u';
    warning('moments_u has been transformed from a row vector to a column vector.');
  end
end

% ENSURE moments_L IS A COLUMN VECTOR
if (size(moments_l,2) > 1)
  if (size(moments_l,1) > 1)
    error('moments_l must be a vector.'); 
  else
    moments_l = moments_l';
    warning('moments_l has been transformed from a row vector to a column vector.');
  end
end

% CHECK DIMENSIONS OF INPUTS TO ENSURE THEY MATCH
J_upper = size(moments_u,1); 
if (size(Omega_u,1) ~= size(Omega_u,2))
  error('Correlation matrix Omega_u must be square.');
elseif (size(Omega_u,1) ~= J_upper)
  error('Number of rows and columns of Omega_u must match number of moments.');
end

% CHECK DIMENSIONS OF INPUTS TO ENSURE THEY MATCH
J_lower = size(moments_l,1); 
if (size(Omega_l,1) ~= size(Omega_l,2))
  error('Correlation matrix Omega_l must be square.');
elseif (size(Omega_u,1) ~= J_upper)
  error('Number of rows and columns of Omega_l must match number of moments.');
end

% Compute analog estimates of lower and upper bounds
bound_estimates = [max(moments_l),min(moments_u)];

% Draw from multivariate normal distribution for computing critical values
sim_draws = normrnd(0,1,max(J_upper,J_lower),R);

% COMPUTE UPPER BOUND OF CONFIDENCE INTERVAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the contact set.
if (strcmp(Selection_Method,'IIV'))
  moment_selection_factor = c * max(se_u);
  contact_set = IIV_moment_selection_upper(moments_u,moment_selection_factor,n);
else
  contact_set = AIS_DiscreteCLR_Upper(moments_u,Omega_u,se_u,n,sim_draws(1:J_upper,:));
end
Omega_contact_set = Omega_u(contact_set,contact_set);

% Determine the critical value.
J = size(contact_set,1);
if (J==1)
  kp = norminv(p,0,1); % Use the standard normal quantile
else
  kp = quantile(max(chol(Omega_contact_set)' * sim_draws(1:J,:),[],1),p);
end
CI_upper_bound = min(moments_u + se_u*kp);
% DONE COMPUTING UPPER BOUND OF CONFIDENCE INTERVAL %%%%%%%%%%%%%%%%%%%%%%%%

% COMPUTE LOWER BOUND OF CONFIDENCE INTERVAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the contact set.
if (strcmp(Selection_Method,'IIV'))
  moment_selection_factor = c * max(se_l);
  contact_set = IIV_moment_selection_lower(moments_l,moment_selection_factor,n);
else
  contact_set = AIS_DiscreteCLR_Lower(moments_l,Omega_l,se_l,n,sim_draws(1:J_lower,:));
end
Omega_contact_set = Omega_l(contact_set,contact_set);

% Determine the critical value.
J = size(contact_set,1);
if (J==1)
  kp = norminv(p,0,1); % Use the standard normal quantile
else
  kp = quantile(max(chol(Omega_contact_set)' * sim_draws(1:J,:),[],1),p);
end
CI_lower_bound = max(moments_l - se_l*kp);
% DONE COMPUTING LOWER BOUND OF CONFIDENCE INTERVAL %%%%%%%%%%%%%%%%%%%%%%%%

confidence_interval = [CI_lower_bound,CI_upper_bound];

end

% Function IIV_moment_selection_lower implements inequality selection for
% lower bound estimates as described on page 666 of NR.
% Function IIV_moment_selection_upper
% implements the analogous procedure for upper bound estimates.
%
% INPUT PARAMETERS:
%   moments: Bound estimates from sample data.
%   c: Multiplicative factor in moment selection. The one described on
%      NR page 666 is 2*max(s_n).
%   n: sample size of the dataset used to compute bound_estimates and
%       correlation_matrix
% OUTPUT PARAMETERS:
%   contact_set: a column vector of indices indicating which entries of
%       bound_estimates are close to the maximum/minimum.

% Computes the contact set for lower bound estimates
function contact_set = IIV_moment_selection_lower(moments,c,n)
  Lstar = max(moments);
  contact_set = find(moments >= (Lstar - c*sqrt(log(n)))*ones(size(moments)));
end

% Computes the contact set for upper bound estimates
function contact_set = IIV_moment_selection_upper(moments,c,n)
  Ustar = min(moments);
  contact_set = find(moments <= (Ustar + c*sqrt(log(n)))*ones(size(moments)));
end

% Function AIS_DiscreteCLR_Lower implements adaptive inequality selection over a
% finite number of inequalities following the method described in
% "Intersection Bounds: Estimation and Inferece," Chernozhukov, Lee, and
% Rosen (Econometrica, 2013).  Function AIS_DiscreteCLR_Upper implements
% the analogous adaptive inequality selection for an upper bound estimate.
%
% INPUT PARAMETERS:
%   bound_estimates: Bound estimates from sample data
%   correlation_matrix: The estimated correlation matrix of the bound estimates.
%   s_n: standard errors of bound estimates
%   n: sample size of the dataset used to compute bound_estimates and
%       correlation_matrix
%   sim_draws: J*R vector matrix of independent N(0,1) draws, where J is
%       number of elements (rows) of bound_estimates
%
% OUTPUT PARAMETERS:
%   contact_set: a column vector of indices indicating which entries of
%       bound_estimates are close to the maximum/minimum, according to the
%       adaptive inequality selection procedure of Chernozhukov, Lee, and
%       Rosen (Econometrica, 2013).
%

% Computes the contact set for lower bound estimates
function contact_set = AIS_DiscreteCLR_Lower(bound_estimates,correlation_matrix,s_n,n,sim_draws)
  p_bar = 1-1/(10*log(n));
  k_bar = quantile(max(chol(correlation_matrix)' * sim_draws,[],1),p_bar); % compute a high-level quantile of the max of the appropriate multivariate normal statistic
  Lstar = max(bound_estimates - k_bar*s_n); % compute the maximum of the precision-corrected estimates
  contact_set = find(bound_estimates >= (Lstar - 2*s_n*k_bar)); % return the indices of bound_estimates that are close to the minimum
end

% Computes contact set for upper bound estimates
function contact_set = AIS_DiscreteCLR_Upper(bound_estimates,correlation_matrix,s_n,n,sim_draws)
  p_bar = 1-1/(10*log(n));
  k_bar = quantile(max(chol(correlation_matrix)' * sim_draws,[],1),p_bar); % compute a high-level quantile of the max of the appropriate multivariate normal statistic
  Ustar = min(bound_estimates + k_bar*s_n); % compute the maximum of the precision-corrected estimates
  contact_set = find(bound_estimates <= (Ustar + 2*s_n*k_bar)); % return the indices of bound_estimates that are close to the minimum
end

